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Abstract:  An  algorithm  -for  selection  of  knot  point  locations  for 

approximation  of  functions  from  large  sets  of  scattered  data  by 
least  squares  Thin  Plate  Splines  is  given.  The  algorithm  is 
based  on  the  idea  that  each  data  point  is  equally  important  in 
defining  the  surface,  which  allows  the  knot  selection  process  to 
be  decoupled  from  the  least  squares.  Properties  of  the  algorithm 
are  investigated,  and  examples  demonstating  it  are  given. 

Results  of  some  least  squares  approx i mat ions are  given  and 
compared  with  other  approximation  methods.  I  ■  ,  * 
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1.0  INTRODUCTION 


The  problem  of  -fitting  a  surface  to  small  sets  of  given  data 
has  been  addressed  in  many  different  ways  and  several  computer 
programs  are  currently  available  which  enable  one  to  deal  with 
the  problem  effectively.  Many  of  the  methods  available  involve  a 
global  1 nterpol ati on  or  approx i mat i on  scheme  and  often  involves 
solving  a  system  of  equations  with  an  equivalent  number  of 
unknowns.  For  very  large  sets  of  data,  the  problem  is 
computationally  intractable.  This  consi derat i on  provides  the 
motivation  behind  the  development  of  a  way  to  pare  the  problem 
down  to  a  more  manageable  size. 

We  wish  to  construct  a  function  F  which  approximately  fits 
the  data  since  we  assume  the  data  collection  is  subject  to 
measurement  error.  We  propose  to  use  approximation  by  least 
squares  Thin  Plate  Splines  (TPS),  where  the  surface  function  is 
constructed  so  as  to  minimize  an  error  function  subject  to 
certain  constraints.  Salving  the  approximation  problem  will  also 
involve  as  many  equations  as  there  are  data  points,  but  the 
number  of  unknowns  will  be  significantly  fewer.  Part  of  the 
appeal  of  TPS  approximation  lies  in  the  fact  that  it  minimizes  a 
certain  linear  functional,  and  involves  a  linear  combination  of 
functions  with  no  greater  complexity  than  the  natural  logarithm 
of  the  distance  function. 

Interpolation  of  scattered  data  by  the  method  of  TPS  was 
developed  from  engineering  considerations  by  Harder  and  Desmarais 
C51 .  It  can  be  thought  of  as  a  two  dimensional  generalization  of 
the  cubic  spline,  which  models  a  thin  beam  under  point  loads 
subject  to  equilibrium  constraints.  The  TPS  function  is  derived 


from  a  differential  equation  which  gives  the  deformation  of  an 
infinite,  thin  plate  under  the  influence  of  point  loads.  A  point 
load  is  applied  at  each  data  point  so  that  the  interpolating 
surface  can  be  constructed  as  a  sum  of  fundamental  solutions  of 
the  TPS  equation. 

In  using  the  least  squares  TPS  approximation  method  to  fit 
the  surface,  a  fewer  number  of  basis  functions  than  the  number  of 
given  data  points  is  employed.  These  basis  functions  are 
centered  at  a  different,  smaller  set  of  points,  which  in  analogy 
with  the  univariate  case,  we  call  the  knots.  Therefore,  the 
problem  at  hand  is  one  of  selecting  the  knot  points,  and  hence 
the  basi s  functions.  This  approach  differs  from  the  use  of 
smoothing  splines,  which  were  introduced  by  Wahba  and 
Wendelberger  C91  in  the  multidimensional  case,  and  called 
Laplacian  Smoothing  Splines  (LSS) .  LSS  minimize  a  certain 
functional  which  is  a  linear  combination  of  a  term  measuring 
fidelity  to  the  data  and  one  measuring  smoothness  of  the  function 
(a  generalization  of  the  usual  thin  plate  spline  functional).  In 
this  case  there  is  still  one  basis  function  for  each  data  point, 
but  the  interpolation  condition  is  relaxed. 

Given  a  'large'  set  of  data  points,  <Xj,yi,fi>,  i  *  1,...,N, 
we  wish  to  find  a  smaller  set  of  knot  points,  (#^,9^),  j  * 
1,...,K,  which  will  'represent'  the  former  reasonably  well.  This 
could  be  accomplished  by  choosing  a  subset  of  the  original  set, 
or  by  some  process  which  produces  a  representat i ve  set.  The 
ultimate  goal  is  to  approximate  the  surface  from  which  the  origi¬ 
nal  data  arose  using  the  representat i ve  set.  Hence,  a  surface 
fit  to  the  large  set  and  one  fit  to  the  representat i ve  set  should 
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Me  use  LINPACK  Cl]  subroutines  to  do  the  actual  calculations. 

Previous  attempts  have  been  made  to  minimize  the  error 
function  by  considering  it  to  be  a  function  of  the  knot  point 
locations  as  well  as  the  coefficients,  wherein  a  total  of  3K 
parameters  are  involved.  As  reported  on  by  Schmidt  [S3,  the 


initial  knot  configuration  was  taken  to  be  of  tensor  product 


form.  The  overall  minimization  process  is  a  large  non-linear 
one,  and  is  complicated  by  possible  coalescense  o-f  knots  as  well 
as  non-unique  solutions  (as  indicated  by  consideration  o-f  one¬ 
dimensional  cases).  Also,  the  objective  function  may  have  many 
local  minima  so  that  avoiding  poor  local  minima  or  searching  for 
better  local  minima  may  be  necessary.  Because  of  these  kinds  of 
problems,  our  goal  is  to  decouple  the  knot  selection  process  from 
the  least  squares  process. 

When  data  are  somewhat  uniformly  distributed,  methods  invol¬ 
ving  tensor  product  cubic  splines  may  be  desirable.  Tensor  pro¬ 
duct  methods  place  knot  locations  on  a  grid,  which  may  not 
reflect  the  actual  disposition  of  the  data  points;  in  fact, 
there  could  be  no  data  nearby.  Furthermore,  even  though  these 
problems  are  surmountable,  they  could  lead  to  nonuniqueness  of 
solutions  and  a  minimum  norm  solution  that  may  not  be 
aesthetically  appealing. 

A  different  point  of  view  is  considered  here  wherein  the 
knot  point  locations  are  predetermined  based  on  two  criteria. 
Specifically,  we  shall  make  assumptions  relating  the  density  of 
data  to  the  dependent  variable  and  mandating  the  importance  of 
each  individual  data  point.  Solution  of  the  overdetermi ned  sys¬ 
tem  of  equations  follows  the  knot  point  selection.  A  summary  of 
the  approach  and  results  will  be  presented.  Examples  are  given 
which  illustrate  rather  well  the  ability  of  the  scheme  to  select 
knot  locations  which  reflect  the  underlying  density  of  the  data. 
Actual  surface  fitting  and  comparison  with  two  other  methods,  the 
Laplacian  Smoothing  Splines  of  Wahba  and  Wendelberger  [93,  and 
the  tensor  product  bicubic  Hermite  method  due  to  Foley  C33  are 
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also  reported  on. 


We  also  point  out  two  related  ideas  which  are  attempts  to 
decrease  the  number  of  basis  functions  to  be  considered.  Each  is 
an  attempt  to  choose  a  subset  of  points  to  be  used  to  construct 
the  approximation.  Schiro  and  Williams  C71  used  an  adaptive 
process  for  subset  selection  and  overlapping  regions  to  construct 
underwater  topographic  maps.  Gazzini ,  diTisi,  and  Lenarduzzi  C21 
gave  an  algorithm  for  selecting  a  subset  of  points  which  were 
important  to  proper  definition  of  the  surface.  Both  of  these 
methods  made  no  assumptions  about  the  density  of  the  data  points 
relative  to  the  behavior  of  the  surface,  and  required 
consideration  of  the  ordinate  values. 

2.0  THE  KNOT  SELECTION  PROCESS 

Given  'a  priori'  flexibility  in  knot  placement,  the  problem 
becomes  the  selection  of  knot  locations,  followed  by  solution  of 
the  system  by  least  squares.  Since  the  selection  of  knot 
locations  is  to  be  decoupled  from  the  solution  of  the  least 
squares  problem,  some  assumptions  must  be  made  in  order  to  devel¬ 
op  an  algorithm  for  the  knot  selection  process.  First,  we 
assume  that  the  independent  variable  data  reflects  something 
about  the  behavior  of  the  dependent  variable.  For  example,  the 
density  of  the  data  points  may  be  dependent  on  the  curvature  of 
the  surface.  Hence,  where  relatively  many  data  points  are  found, 
the  function  is  assumed  to  be  changing  behavior  rapidly,  whereas 
a  low  density  of  data  indicates  slowly  changing  behavior. 

Although  this  assumption  is  not  universally  satisfied  in 
practice,  it  does  not  seem  to  be  an  unreasonable  one. 


The  second  assumption  is  that  each  data  point  is  equally 
important  in  defining  the  underlying  surface.  Therefore  the 
number  of  data  paints  represented  by  each  knot  should  be  the  same 
or  nearly  the  same.  This  leads  to  equal  representation '  of  the 
data  points  by  the  knot  points  where  each  data  point  is  close' 
to  a  knot  point.  A  key  advantage  achieved  in  pursuing  this 
approach  is  the  existence  of  a  natural  heuristic  for  moving  the 
knots  around  the  plane  in  searching  for  a  good  knot 
configuration.  This  point  will  be  elaborated  on  later  in  the 
paper . 

Our  knot  selection  algorithm  is  based  on  these  last  two 
assumptions.  First,  we  wish  to  minimize  the  sum  of  the  distances 
squared  from  each  data  point  to  the  nearest  knot  point;  that  is, 
minimize  the  'global'  value, 

,  N 

GN^  =  £  min  t (x . -k  )2  +  <y 

1=1  j  1  J  J 

This  function  is  continuous  and  piecewise  quadratic.^  The 
expression  leads  naturally  to  a  'default'  Dirichlet  Tesselation, 
a  partitioning  of  the  plane  with  respect  to  the  knot  points  (see 
Figure  1).  Thus,  we  say  each  data  point  belongs  to  some  knot 
point  according  to  the  Dirichlet  tile  in  which  it  lies.  Data 
paints  on  any  of  the  tile  boundaries  (ties)  must  be  resolved  by  a 
determination  of  which  tile  they  belong  to  or  some  sharing 
mechanism.  Our  initial  guess  at  the  knot  point  configuration  was 
taken  to  be  quasi -gri dded. 

2 

Differentiation  of  GN  with  respect  to  the  and  show 
that  at  the  minimum,  each  knot  point  will  occupy  the  centroid 
with  respect  to  the  data  points  inside  that  tile.  Given  the 


initial  configuration  of  knot  points  with  its  Dirichlet 


Tesselation,  the  following  algorithm  for  iteration  to  a  local 
2 

minimum  GN“  value  is  employed: 

<a)  compute  the  centroid  of  each  tile  with  respect  to  the 
data  paints  contained  within  each  tile; 

<b)  move  the  knots  to  the  corresponding  centroids,  which 
results  in  a  new  Dirichlet  Tesselation  and  a  new  set  of  knot 
point  -  data  point  associations;  this  is  the  configuration  for 
the  next  iteration. 

(c)  quit  when  two  successive  iterations  yield  the  same  knot 
locations,  which  means  that  a  local  minimum  value  of  Gn  has  been 
found. 

This  algorithm  was  formulated  in  discussions  at  the  Istituto  per 

le  Applicazioni  della  Matematica  e  del  1  ' Inf ormat i ca  in  1983  C103, 

after  the  problem  was  posed  by  G.  Nielson  and  R.  Franke. 

2 

We  note  that  the  value  of  GN  will  necessarily  decrease  as 
the  iterations  continue,  until  two  successive  iterations  yield 
the  same  conf iguration;  this  will  be  proven  below.  In  the  case 
where  no  data  points  lie  in  a  tile  for  some  knot  point,  the  knot 
point  is  moved  to  the  nearest  data  point.  This  mechanism  avoids 
knots  without  data  points.  Futhermore,  if  a  data  point  lies  on  a 
tile  boundary,  it  is  assigned  to  the  knot  with  the  smallest 
subscript  (amongst  the  appropriate  choices  of  knot  points). 
Employment  of  a  different  criterion  for  the  resolution  of  ties 
may  yield  different  results.  We  note  that  knots  cannot  coalesce. 

The  following  theorem  is  pertinent  to  this  algorithm. 

2 

Theorem:  The  function  GN  decreases  with  each  iteration  which 

involves  movement  of  a  knot  point. 

2 

Proof:  Write  GN  in  the  more  convenient  form 

K 

GN2  =  Z  Z  C  (x  .  — {f  . ) 2  +  (y.-y.)2]  (1) 

j=l  i« I ^  1  J  1  J 
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where  I 


=  Cl:  <x  i  » V  j  )  in  the  tile  -for  >>.  In  (1),  the 

interior  sum  is  the  sum  of  the  distances  squared  from  the  data 
points  in  a  tile  to  the  knot  point  in  that  tile,  and  the  exterior 
sum  is  over  all  K  of  the  tiles.  Let  a  prime  denote  the  new  knot 
points  and  index  sets.  This  form  leads  to  the  expressions 

=  <  Z  x./n.,  Z  y./n.), 

J  i‘  I  J  ifi I  .  J 

J  J 

where  n^  is  the  number  of  indices  in  the  set  1^.  The  new  knot 
points  will  lead  to  a  new  tesselation,  followed  by  the  new  index 
sets  IJ.  Then  the  expression  (1)  is  greater  than  or  equal  to 

K  -  2 

Z  E  C  (x. -*'.)*  ♦  (y.-y.)  1  (2) 

j  =  l  id.  1  J  1  J 

because  the  new  knot  point  locations  minimize  the  contribution  of 
the  interior  sums.  This  expression  (2),  in  turn,  is  greater  than 
or  equal  to 

K  2  2 

z  Z  C  <x  -#' )  +  (y .  — y '.  >  j  (3) 

j=l  iclj  J  1  J 

since  an  index  i  moves  to  another  set  only  in  the  case  wherein 

the  corresponding  data  point  is  now  closer  to  a  different  knot 

2 

point,  thus  decreasing  its  contribution  to  the  global  GN  value. 

2 

Finding  a  local  minimum  of  GN  is  well-served  by  this  algo¬ 
rithm;  however,  as  seen  in  a  one  dimensional  example  C61,  the 
2 

function  GN  is  rife  with  local  minima,  and  the  local  minimum 
value  found  depends  on  the  initial  configuration  of  knots  used. 

We  can  draw  similar  conclusions  for  the  multidimensional  case 
based  on  the  one  dimensional  analogy. 

The  process  of  locating  each  knot  occurs  in  two  distinct 
steps:  first,  each  data  point  is  assigned  to  the  closest  knot 


and  second,  a  determination  is  made  within  each  tile  as  to  the 

location  of  its  centroid.  As  a  direct  result  of  the  centroid 

2 

requirement,  the  GN  function  will  stabilize  at  the  local  minimum 
value  corresponding  to  the  particular  quadratic  piece  on  which 
the  knot  is  found. 

The  local  minimum  value  will  frequently  occur  out  of  the 
domain  of  the  correspondi ng  quadratic  piece.  This  leads  to  a 
'cascading'  phenomenon  which  continues  until  a  local  minimum 
value  is  attained.  However,  the  global  minimum  value  will  not 
necessarily  be  attained  since  the  cascading  will  stop  as  soon  as 
a  local  minimum  value  is  found  within  the  domain  of  each 
quadratic  piece. 

This  inconsistent  performance  of  the  algorithm  in  finding 

2 

the  global  minimum  value  of  GN  leads  to  consi derati on  of  a 
somewhat  different  criterion  for  locating  a  good  configuration  of 

knot  points.  We  wish  to  exploit  the  second  assumption  specified 

2 

earlier,  while  taking  advantage  of  the  minimization  of  the  GN 
function.  Since  each  data  point  is  assumed  to  be  equally  impor — 
tant,  the  Dirichlet  tile  for  each  knot  should  contain  about  the 
same  number  of  data  points.  Thus,  we  wish  to  minimize  the  sum  of 
the  squares  of  the  differences  between  the  number  of  knots  in 
each  tile  and  the  average  number  of  data  points  that  should 
belong  to  each  tile;  that  is,  minimize  the  quantity 

K  2 
D  =  Z  <n  -  N/K) 

j  =  l  J 

The  new  algorithm  for  determining  knot  locations  is  based  on  the 
minimization  of  D,  subject  to  the  constraint  that  each  knot  be 
located  at  the  centroid  of  its  tile. 
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k 


This  new  optimization  leads  to  a  natural  heuristic  tor 
moving  knots  from  a  stable  configuration  to  a  possibly  better 
con-figuration.  We  call  the  current  configuration  of  knots  a 
'base  configuration,  and  iterate  through  the  algorithm  as 
foil ows: 

(a)  generate  a  new  guess  for  the  knot  locations  by  moving 
the  knot (s)  with  the  smallest  number  of  data  points  in  their 
tile(s)  toward  the  knot(s)  with  the  largest  number  of  knot (s)  in 
their  tile(s);  the  distance  moved  is  initally  a  large  fraction  of 
the  total  distance  between  the  knots. 

(b)  iterate  to  a  stable  configuration  using  the  first 
algorithm,  compute  the  values  of  GN^  and  D,  and  compare  D  to  the 
smallest  value  achieved  to  date,  as  represented  by  that  of  the 
base  conf i gurat i on ; 

(c>  repeat  the  process  above  when  a  smaller  value  of  0  is  ob¬ 
tained,  with  the  present  configuration  as  the  base  configuration; 

(d>  when  a  smaller  value  of  D  is  not  found,  take  a  shorter 
step  in  the  movement  of  the  knot(s)  and  repeat  the  process  above; 

(e)  continue  with  smaller  and  smaller  steps  until  a  smaller 
value  of  D  is  found  (or  an  equal  value  of  D  with  a  smaller  Glvr 
value)  or  until  the  knot  locations  return  to  the  base 

conf 1 gurat ion; 

(f)  perform  the  search  in  the  symmetrical  way  when  the  base 
configuration  is  returned  to;  that  is,  move  the  knot(s)  with  the 
largest  number  of  data  points  in  their  tiles  toward  the  knot(s) 
with  the  smallest  number  of  points  in  their  tile; 

<g)  quit  when  no  smaller  value  of  D  is  found. 

The  movement  of  the  knots  is  justified  by  the  rationale  that 

a  more  equitable  distribution  of  data  points  can  be  found  by 

moving  the  tile  boundaries  across  data  points.  Note  that  the 

2 

algorithm  for  computing  a  local  minimum  of  the  GN  function  value 
is  embedded  in  this  new  algorithm. 


3.0  RESULTS  AND  EXAMPLES 

Using  our  algorithm  for  the  a  priori  selection  of  the  knot 
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point  locations,  experiments  were  conducted  to  test  the  scheme 


using  different  sets  of  test  data.  This  was  followed  by 
verification  of  the  scheme  on  a  large  set  of  real  data.  Results 
from  two  sets  of  the  test  data  are  presented  here:  one 
consisting  of  200  data  points  called  Cliff',  and  one  consisting 
of  500  data  points  called  'Humps  and  Dips'.  Both  sets  of  data 
were  generated  using  known  functions  (see  Franke  C4D)  in  a  way 
that  forced  the  density  of  points  to  be  proportional  to  the 
curvature  of  the  sampled  function.  Figures  2-5  show  these  two 
test  data  sets  graphically,  and  illustrate  the  optimized  knot 
point  configurations  found  using  the  least  squares  algorithm. 
Figure  6  depicts  actual  hydrographic  data  collected  in  Monterey 
Bay,  with  greatly  varying  density.  Figure  7  shows  the  results  of 
applying  the  algorithm  with  100  knots,  and  is  particularly 
interesting  because  the  density  of  the  data  is  faithfully 
replicated  by  the  knots.  Ule  note  that  the  assumption  regarding 
the  density  of  the  data  points  being  indicative  of  the  behavior 
of  the  dependent  variable  is  not  actually  satisfied  in  this  case, 
due  to  the  source  of  the  data.  Nonetheless,  these  results 
demonstrate  the  ability  of  the  algorithm  to  produce  representa¬ 
tive  sets  of  knots. 

Me  also  investigated  how  closely  the  constructed  surface  F 
and  the  'true*  surface  resemble  one  another.  This  comparison  is 
made  in  the  context  of  the  root-mean-squared  error  (RMS)  of  the 
residuals  (at  the  data  points)  and  on  a  rectangular  grid  of 
locations  in  the  plane.  The  two  data  sets  constructed  above,  and 
one  other  which  was  generated  in  a  similar  manner,  were  used  to 
compare  RMS  errors  for  the  least  squares  algorithm  developed  here 

1 1 


with  the  method  of  Wahba  and  Wendelberger  (Laplacian  Smoothing 
Splines),  and  the  method  of  Foley  (Bicubic  Hermite  Tensor  Product 
Spl ines) . 

The  dependent  variable  values  of  the  experimental  data  sets 
were  generated  in  two  ways:  1>  using  a  known  function,  and  2) 

contaminating  the  fcnown  function  by  the  injection  of  independent, 
identically  distributed  normal  random  errors  with  a  specified 
composite  standard  deviation  of  0.05.  The  actual  standard 
deviation  was  about  0.0485.  In  the  first  case,  we  would  expect 
the  RMS  error  on  the  data  points  and  on  the  grid  to  be  about  the 
same,  and  to  decrease  as  the  number  of  knot  points  used  to 
represent  the  data  is  increased. 

In  the  contaminated  case,  the  dependent  variable  at  each  data 
point  is  the  sum  of  the  unknown  underlying  function  value  and  the 
error  function  value  so  that  the  difference  between  the 
constructed  surface  and  the  'true'  surface  is  mainly  attributable 
to  the  presence  of  error  in  the  data.  In  the  best  situation,  we 
expect  the  RMS  error  in  the  residuals  to  match  the  composite 
standard  deviation  of  the  random  error  injected  to  obtain  the 
contaminated  data.  At  the  grid  points,  we  expect  the  RMS  error 
to  be  smaller  than  the  composite  standard  deviation,  since  the 
grid  sample  is  larger  (33*33)  and  the  errors  are  distributed  more 
evenly  throughout  the  entire  region  of  interest. 

Some  observations  can  be  made  regarding  Tables  1-3.  The 
general  trend  of  the  RMS  error  on  both  the  data  points  and  the 
grid  is  to  diminish  as  the  number  of  knot  points  is  increased. 

As  expected  with  the  exact  data,  the  RMS  error  of  the  residuals 
and  the  RMS  error  on  the  grid  are  roughly  equivalent.  For  the 


contaminated  data,  the  RMS  error  of  the  residuals  roughly  matches 
the  composite  standard  deviation  o-f  the  data,  and  the  RMS  error 
on  the  grid  is  smaller  than  the  RMS  error  o-f  the  residuals,  as 
expected.  In  Table  1,  the  errors  over  the  grid  increase  as  the 
number  of  knot  points  is  increased,  and  the  RMS  value  of  the 
residuals  becomes  less  than  the  composite  standard  deviation  of 
the  injected  errors.  In  this  case,  under smoothi ng  has  occurred, 
and  the  surface  is  "drawing"  toward  the  error. 

In  comparing  least  squares  TPS  to  the  smoothing  spline 
method  in  the  exact  data  case,  we  note  that  the  smoothing  spline 
method  yields  a  residual  RMS  error  of  0.  This  could  be  expected, 
since  there  is  no  error  in  the  data  and  the  spline  of 
interpolation  is  chosen.  On  the  grid,  the  RMS  error  is  small. 
When  the  data  is  not  contaminated,  the  RMS  error  of  the  least 
squares  algorithm  only  begins  to  become  as  small  as  that  of  the 
smoothing  splines  method  when  the  number  of  knots  used  becomes 
large.  We  also  note  that  in  the  500  data  point  set,  no 
comparison  is  made  since  a  potential  limit  for  computing 
smoothing  splines  is  200-300  data  points. 

Foley's  method  for  the  contaminated  case  gives  errors  nearly 
equal  to  the  composite  standard  deviation  of  the  errors  injected 
into  the  data.  However,  on  the  grid,  the  least  squares  method 
does  better,  an  indication  that  smoothing  is  occurring,  as 
expected.  We  also  note  that  an  increase  in  the  number  of  grid 
points  does  not  significantly  improve  the  RMS  error  in  Foley's 
method,  even  though  an  increase  in  the  number  of  knots  in  the 
least  squares  method  usually  yields  improved  results.  We  used 
the  default  local  approximations  in  Foley's  method,  and  we  note 
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that  performance  of  the  method  may  be  improved  by  using  lower 
degree  local  approximations  to  estimate  the  grid  values  to  be 
used . 

Finally,  we  note  that  the  search  for  a  best  knot  configura¬ 
tion  can  turn  out  to  be  rather  expensive.  For  a  large  number  of 
data  points  with  a  moderately  large  number  of  knot  points,  the 
computat i anal  effort  could  be  excessive,  although  we  are 
investigating  ways  of  speeding  up  the  algorithm.  Furthermore,  as 
we  noted  earlier,  the  end  results  are  dependent  on  the  initial 
guess,  although  they  generally  look  quite  good  for  any  reasonable 
initial  guess. 
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Figure  Is  A  Dirichlet  Teeselation  with  5  Tiles.  It  is 
constructed  by  connecting  the  perpendi cul ar  bisectors  of  the 
lines  joining  each  of  the  knot  points. 


Figures  4  and  5:  The  Humps  and  Dips'  data  set.  Note  how 

clumps  of  data  appear  in  three  portions  of  the  plane,  indicating 
that  the  underlying  surface  is  undergoing  change.  A  set  of  jU 
I  not  points  was  used  to  represent  the  delta. 
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